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Abstract. New numerical models of line-driven stellar winds of late O stars are presented. Statistical equilibrium (NUTE) 
equations of the most abundant elements are solved. Properly obtained occupation numbers are used to calculate consistent 
radiative force and radiative heating terms. Wind density, velocity and temperature are calculated as a solution of model hydro- 
dynamical equations. Contrary to other published models we account for a multicomponent wind nature and do not simplify 
the calculation of the radiative force (e.g. using force multipliers). We discuss the convergence behaviour of our models. The 
ability of our models to predict correct values of mass-loss rates and terminal velocities of selected late O stars (mainly gi- 
ants and supergiants) is demonstrated. The systematic difference between predicted and observed terminal velocities reported 
in the literature has been removed. Moreover, we found good agreement between the theoretical wind momentum-luminosity 
relationship and the observed one for Cyg 0B2 supergiants. 
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1. Introduction 

NLTE line-driven stellar wind models have been calculated by 
various authors (cf. de Koter et al. 1993' Hillier & Miller 1998 
Vink et al. ,.2000. Pauldrach et al. 2001 Griifener et al. 2002 1. 
These sophisticated models cover a wide range of OB and WR 
star parameters, and it may seem that there is no need for any 
additional NLTE wind models. However, we believe this is not 
the case. Firstly, these models employ different simplifying as- 
sumptions which cannot be necessarily fulfilled in all model 
stars. For example, wind velocity is often calculated using ei- 
ther a simplified /3- velocity law (Vink et al. 120011 Grafener et 
al. I2002> or using simplified line-force multipliers (Pauldrach 
et al. l200n Kudritzki 2002 1 of Castor et al. (1975 J and Puis 
et al. (l2000t . Some of these models (e.g., Vink et al. ISOOH 
do not allow independent calculation of observed terminal ve- 
locities. Moreover, temperature stratification is often neglected 
(cf. Kudritzki l2002l assuming that the radiative force is not 
significantly influenced by temperature (see Pauldrach '1987i 
or calculated using a simplified expression (Pauldrach et al. 
119941 Vink et al. I1999> . Although the calculation of correct 
temperature stratification is included in modern wind codes (cf. 
Pauldrach et al. ^01 1, the latest published grid of NLTE cal- 
culations of wind temperature dates back to the paper of Drew 
J1989> . who used a relatively constrained set of statistical equi- 
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librium equations. Despite the abovementioned enormous ef- 
fort in wind model calculations, many unsolved problems re- 
main which may be connected either with the simplifying as- 
sumptions mentioned above or with some others, e.g. neglect 
of wind non-stationarity (cf. Runacres & Owocki '20021, wind 
3D structure (Dessart & Owocki 2003 1, or wind magnetic fields 
(ud-Doula & Owocki |2002> . Among these unsolved problems 
is the long-standing momentum problem of WR stars (for the 
latest progress see Grafener et al. 2002 ) or problems with fitting 
of X-ray line profiles (Kramer et al. l2003l l. 

Many of the simplifying assumptions employed in mod- 
ern wind codes are probably acceptable for high density 
radiatively-driven stellar winds of O stars, but they may break 
down for low density stellar winds. Whereas the bulk wind 
material consists of hydrogen and helium, line-driven stellar 
winds are accelerated mainly by the absorption of radiation in 
the resonance lines of trace elements like C, N, O, or Fe. For 
low-density stellar winds this multicomponent nature may be 
important for the wind structure itself, influencing both tem- 
perature stratification (Springmann & Pauldrach T99T Krticka 
& Kubat 2001 hereafter KKII) and wind dynamics (Babel 
[19951 Porter & Skouz a [T999I KKII, Owocki & Puis |2002| 
Krticka & Kubat 2002). Clearly, such a complicated situation 
cannot be adequately described by models with a temperature- 
independent radiative force, or with a radiative force artificially 
depending on wind temperature via the thermal velocity as has 
been done, e.g., by KKII. Note that the artificial dependence 
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on the thermal velocity may be removed using an alternative 
description of the radiative force by Gayley J1995> . 

The situation is especially appealing now, when wind mod- 
els suitably describing the atmospheres of first generation stars 
have been calculated (Kudritzki 2002 1, because some of these 
models are also influenced by multicomponent effects (Krticka 
et al. 120031 see also Kudritzki 2002.) . Moreover, predicted 
mass-loss rates of low-metallicity winds of some SMC stars 
are more than one magnitude higher than the observed ones 
(Bouret et al. .2003i) . These discrepancies may be connected 
with problems with source function gradients (Owocki & Puis 
[T999I . 

Although satisfactory determination of mass-loss rates and 
terminal velocities of O stars can be performed from observa- 
tions in different wavelength regions, the situation in the B star 
domain is not so clear. Even for massive B stars there is a sig- 
nificant discrepancy between observations in different wave- 
length regions (Vink et al. 2000 1. Even worse, for most main- 
sequence B stars there is no trustworthy determination of either 
mass-loss rates or terminal velocities, and, moreover, it is not 
clear which B stars posses a stellar wind (see Babel ( 1996) for 
the discussion of this problem). This is a challenge for sophis- 
ticated wind models to help supplement missing observational 
data. There exist several subclasses of B stars whose peculiar 
behaviour may be connected with their stellar wind. This is 
epecially the case of Be stars, because the plausible explana- 
tion of the Be phenomenon is still missing (excluding some 
binaries, for which a binary hypothesis of Kfiz & Harmanec 
J1975> offers a trustworthy explanation). 

Last but not least, there is an ambitious idea of the 
wind momentum-luminosity relationship (see Kudritzki & Puis 
<2000> and references therein) which may be used as an inde- 
pendent reliable method for determination of stellar distances. 
Thus, several types of detailed wind models are necessary to 
verify all of the assumptions of line-driven wind models. 

Motivated by these considerations we offer a new indepen- 
dent type of NLTE line-driven wind model to obtain correct 
values of radiative force, mass-loss rates and terminal veloci- 
ties in the case of low-density stellar winds. 

2. Model assumptions 

Hydrodynamic model equations solved in this paper are essen- 
tially the same as in KKII. An interested reader will find more 
details about the solution of multicomponent hydrodynamical 
equations in KKII and Krticka (2003). Models presented here 
differ in the method of calculation of the radiative force and 
the radiative heating/cooling term, and these new models as- 
sume slightly different boundary velocities. In the following 
sections we summarize our model equations. The basic model 
assumptions are the following: 

- We assume a stationary spherically-symmetric flow. 

- We solve the continuity equation, the momentum equation 
and the energy equation for each component of the flow, 
namely for absorbing ions, nonabsorbing ions (hydrogen 
and helium), and electrons. Nonabsorbing ions are ions for 
which the effect of a radiative force may be neglected (e.g.. 



hydrogen and helium if Population I/II stars are consid- 
ered). 

- Level occupation numbers of model ions are obtained from 
NLTE rate equations with an inclusion of a superlevel con- 
cept (Anderson 1989 1. 

- Radiation transfer is split into two parts. Line radiative 
transfer is solved in a Sobolev approximation (Rybicki 
& Hummer 1 1978> neglecting continuum opacity sources, 
line overlaps, and multiple scattering. Continuum radiative 
transfer is formally solved using the Feautrier method with 
an inclusion of all free-free and bound-free transitions of 
model ions, however neglecting Une transitions. Thus, we 
neglect line-blocking. 

- For the solution of NLTE rate equations together with the 
radiative transfer equation we apply the Newton-Raphson 
method for the radiative transfer in lines and lambda iter- 
ations in continuum. The simplification of a lambda itera- 
tion may, of course, fail for optically thick continua, like 
the Hell continuum (however, see Sects. |231 and l3l for a 
discussion of this issue). 

- Derived occupation numbers are used to calculate radiative 
force in the hydrodynamic equations. 

- The radiative cooling/heating term is obtained by adding 
the contributions from all free-free, bound-free and bound- 
bound transitions of model ions. For the calculation of 
this term we use the thermal balance of electrons method 
(Kubat et al.TWl. 

- The description of a transition region between a quasi-static 
stellar atmosphere and supersonic stellar wind is simpli- 
fied. This region is important for the continuum formation 
(Pauldrach_1987 1. This subsonic region was not included 
in our previous models since we did not solve the radia- 
tive transfer equation (see KKII). To obtain correct ionizing 
flux we included this subsonic region in our models, i.e. the 
boundary velocity at the stellar surface is set to some low 
value well below the sonic speed. The particular value of 
the base velocity is selected in such a way that lowering 
of this value by a factor of 10 does not change the model 
structure (i.e. mass-loss rate and terminal velocity) by more 
than 10%. For most of the wind models this means that the 
boundary velocity is of the order of 10"^ of the hydrogen 
thermal speed. 

- The calculation of models below and above the critical 
point is split in two parts, thus we neglect the influence of 
the radiation emitted in the region above the critical point 
on the region below the critical point. However, the influ- 
ence of the subcritical parts on the supercritical ones is 
taken into account coiTectly. 

- We neglect Gayley-Owocki heating (Gayley & Owocki 
I1994> which is not important for stars considered in this 
paper (see KKII). 

- We also neglect the X-ray radiation which may influence 
the ionization balance. 



These assumptions are discussed in more detail in the Sect-ls] 
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2. 1. Equations of statistical equilibrium 

Our treatment of statistical equilibrium equations is very simi- 



lar to that of Pauldrach J1987> . We solve rate equations in the 
form of (Mihalas fl978t 



0, 



(1) 



where Ni is the number density of ions in state i relative to 
the total number density of given atom, and Py , Pji are rates 
of aU processes that transfer an atom from a given state i 
to state j and back, respectively. They are sums of radiative 
and collisional rates (both excitation/de-excitation and ioniza- 
tion/recombination). The system of rate equations Q is closed 
by the equation of normalisation of relative number densities 
for each model atom (i.e. abundance equation) 



(2) 



This equation can be used instead of any equation from the 
set ([0. However, it is computationally convenient to replace 
the rate equation for the level with the highest relative number 
density. 

We use the relative number density Ni instead of the com- 
monly used number density to be consistent with the hydro- 
dynamical part of the code, because relative quantities are ad- 
vantageous from the numerical point of view. Note that relative 
number densities were used for the solution of the equations of 
statistical equilibrium by Lucy feOOlj . 

2.1 .1 . Model atoms and rates 

Radiative excitation rates Rij are given by 



(3) 



where Jy = <f)ij{i')Ji,di' is profile-weighed line intensity 
and 

_j7r2e^ 

nVijlTleC 

where Vij is line frequency, is corresponding oscillator 
strength, e is elementary charge and TOo is electron mass. 
Similarly, the radiative de-excitation rates Rji are 



NjRji = Nj{Aj,+Bj,J,j), 
where the remaining Einstein coefficients are given by 
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and gi, gj are statistical weights of involved levels. 

Radiative ionization and recombination rates are given by 
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Table 1. Atoms and their ionization stages included in the 
NLTE calculations. In this table, level means either an indi- 
vidual level or a set of levels merged into a superlevel. 
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where Ui^u is the corresponding photoionization cross-section, 
Ji, is the mean continuum intensity (line absorption spectrum 
is neglected in calculations of bound-free rates), Vi is the ion- 
ization frequency of level i. To is the electron temperature, and 
an asterisk denotes LTE values. 

Upward (C^) and downward {Cji) collisional rates (for 
both excitation and ionization) are calculated using the expres- 
sions 

N.C^j = N,n,n,j{T,), (8a) 
N,C,,^N,n,(^Eyj n,j{T,), (8b) 

where fi^ (To) is the integrated collisional cross-section (col- 
lision strength) and nc is the electron number density. Finally, 
rate coefficients for both excitation and ionization are sums 
of radiative and collisional terms 



Ri 



(9) 



We solve NLTE rate equations Q for each of the elements 
listed in Table Q separately. This allows us (in addition to a 
fixed temperature) to fix the number of free electrons during 
the solution of the equations of statistical equilibrium. This step 
greatly reduces the computational needs, since instead of one 
large rate matrix we have to deal with a moderate number of 
smaller matrices for the elements considered. 

The model atoms are predominantly taken from TLUSTY 
models (Hubenv .l988t Hubeny & Lanz ,1992. Hubeny & Lanz 
119951 Lanz & Hubenv l2003t . The set of model atoms is slightly 
modified for the case of hot star wind modelling. For this pur- 
pose we used data from the Opacity Project (Seaton il987. Luo 
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& Pradhan [T989l Sawey & BemngtonTOQ^l Seaton et al. [T992l 
Butler et al. [T993l Nahar & Pradhan 1993 1 and from the Iron 
Project (Hummer et al. 1993 Bautista 1996 Nahar & Pradhan 
[TW6i Zhang . 1996 Bautista & Pradhan . 1997. Zhang & Pradhan 
[19971 Chen & Pradhan 1999 1. 

We use a similar strategy for the preparation of model 
atoms that enter the calculations as was akeady applied in 
TLUSTY. With the exception of Fc and Ni, levels with low 
energies (relative to the ground level of a given ionization 
stage) are treated individually (neglecting the hyperfine split- 
ting), whereas higher levels with similar properties are merged 
into superlevels. For Fe and Ni all levels (except the ground 
ones) are merged into superlevels. More details on the con- 
struction of model atoms and atomic data used are given in 
AppendixlAl 

Procedures for the calculation of fiy (Tc) are taken from 
the TLUSTY code. The partition function approximations are 
mostly from Smith & Dworetsky (1988 1 and Smith (1992i. 
Moreover, the partition function approximations taken from the 
TLUSTY code are applied for C, N, O, Fe VI, and Fe VII. 

2.1.2. Superlevels 

Compared to the pioneering work of Anderson (1989), we use a 
simplified merging of individual levels into superlevels. Given 
the sub levels / with energies En, relative number densities Nu, 
and statistical weights gu, the total population of superlevel is 
given by 

iV^=^iV^^, (10) 

I 



and the energy of superlevel is 



(11) 



Thus, we effectively neglect the exp{—Eii/kTc) term. 
Similarly, the photoionization cross section of superlevel a^^i, 
is calculated by 

J2i gaOiiUw 



Y.I 9^l 



(12) 



and the relative number density of a sublevel of a given super- 
level is approximated as 



Na = 



(13) 



Numerical tests showed that the neglect of the exp{—Eii /kT^) 
term does not have a significant effect on the calculated radia- 
tive force. 

2.2. Radiative transfer 

Radiative flux at the base of the stellar wind He is taken 
from H-He spherically symmetric NLTE model atmospheres 
of Kubat (2003 1. The line transfer equation and continuum 
transfer equation are solved separately. 



2.2.1 . Line radiative transfer 

The averaged line intensity is calculated using the Sobolev ap- 
proximation (Rybicki & Hummer 1978 1 



4- = (l-/3)5,, +/3c/c, 



(14) 



where escape probability functions (i and (3^ are given by inte- 
grals 



dfi- 



1 - e-^f^ 



(15) 
(16) 



= (1 - i?2/r2)^/^ and 4 = 4iJc. The Sobolev optica l 
depth is given by (Castor [T974l Rybicki & Hummer fTWSb 
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where Vri is the radial velocity of absorbing ions, rii, nj are the 
number densities of individual states. 



tli = iViZatom"H 



(18) 



(similarly for rij), where Zatom — "-atom/^-H is the number 
density of a given atom relative to the hydrogen number density 

Pp 



= , (19) 

where pp is the density of a passive component (hydrogen and 
helium). The variable cr was introduced by Castor J1974> . 



dlnw,.; ^ 
dlnr 



(20) 



A suitable method for the calculation of escape probability 
functions [3 and (3c was described by KKII. Finally, the line 
source function is 



(21) 



2.2.2. Continuum radiative transfer 

For the calculation of the continuum intensity we use a 
procedure for the solution of the transfer equation based 
on the Feautrier method in the spherical coordinates (for a 
description, see Mihalas & Hummer [T974I or Kubat [1993}. 
However, for the solution of the obtained tridiagonal system 
of equations we use the numerical package LAPACK 
l |http : //www, cs ■ Colorado . edu/ ~lapack| 
Anderson et al. I1999> instead of the classical Gaussian 
scheme. All bound-free transitions of explicit ions, free-free 
transitions of hydrogen and helium and light scattering due to 
free electrons contribute to the source function. The radiative 
transfer equation is solved only above the hydrogen ionization 
edge. For frequencies lower than the hydrogen ionization limit 
we use an optically thin approximation Jvir) — 'iW{r)Hc, 
where W{r) — ^ (1 — /x,) is the dilution factor. 
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2.3. Solution of NLTE equations 

For the solution of the NLTE rate equations Q, Q together 
with the radiative transfer equation we use a combination of 
lambda iterations and of a complete linearization. For the con- 
tinuum transfer we use simple lambda iterations. Although the 
stellar wind is optically thick near to the star below the He II 
ionization limit, these iterations converge relatively well (see 
also Sect. |3}. In lines, we use complete linearization of the 
radiative transfer equation in a simple Sobolev approximation 
(neglecting continuum radiation, see Eq. il4l ). The line inten- 
sity Jij depends only on the level populations of the given atom 
for this case. Thus, the statistical equilibrium equations for dif- 
ferent atoms are independent (for a. fixed hydrodynamical struc- 
ture). Hence, the following scheme may be obtained from sta- 
tistical equilibrium equations Q, (|2} for the calculation of the 
correction 6Ni of relative number densities in the given itera- 
tion step 



E 



BP 



ON,, 



p 



-E 



ON, 



ON. 



= 0, (22) 



where the terms dPij/dNj are calculated using Eqs. il4l - 
i21\ . For details see Appendix |b1 NLTE equations i22\ are 
solved for each element separately. Note that Eq. i22\ is an 
application of a Newton-Raphson method to statistical equi- 
librium equations Q. For the solution of the system of NLTE 
equations ( I22t we use the numerical package LAPACK. Note 
that we can again take full advantage of LAPACK subrou- 
tines designed for general band matrices because our NLTE 
rate equation matrices have a band structure. 

2.4. Derivatives of occupation numbers 

The hydrodynamic equations are solved using a Newton- 
Raphson method. To obtain a well converged model we have 
to know the explicit values of derivatives of the radiative force 
and of the radiative heating/cooling term with respect to fluid 
variables (i.e. densities, temperatures, etc.). However, because 
the radiative force and the radiative heating/cooling term de- 
pend on level populations, we have to calculate derivatives of 
level populations with respect to fluid variables. These can be 
calculated in the following way. For clarity, we rewrite NLTE 
rate equations Q in a symbolic form 

J2Pd^^MN,h),h)N,- 

-iV,^P„ (Ar, J,(Ar,/i),/i) = 0, (23) 

where we denoted all explicit dependencies of rate equations. 
The elements of vector N are level populations Ni, the ele- 
ments of vector h are fluid variables h — (rto, Pp, T^, Vri), and 
the mean continuum intensity J,y depend on h and TV via the 
radiative transfer equation. Clearly, it is not feasible to obtain 
correct exphcit values of dNi / dh (where h denotes elements 



of vector h), since the level populations depend on the solu- 
tion of radiative transfer equation in continuum. However, it is 
possible to obtain approximate values of these derivatives ne- 
glecting the explicit dependence of rates Pij on the continuum 
intensity. Deriving Eq. ( I23t with respect to h and neglecting the 
dependence on Jiy{N, h) we obtain a system of equations for 
dNjdh 



E 
E 



iV, 

ON, dh ^ 



P. 



dh 



dNj dh 



dN,, 



'N,, + P^ 



dh 

m 

dh 



0. (24) 



All calculated derivatives dNi/dh are inserted into the 
Newton-Raphson matrix of hydrodynamic equations (corre- 
sponding to the radiative force, radiative heating/cooling term 



and ionization charge, see KKII and Krticka 1200 H . 

Note that the system of equations for derivatives J24> has 
the same form as the system of NLTE rate equations i22i . 
however with a different right-hand side. Thus, because the 
LAPACK package uses an LU decomposition for the solution 
of the system of linear equations, the values of dNi/dh can 
be calculated using data from the previous solution of NLTE 
equations ( I22t . 

This method is especially useful in the case when the in- 
tensity of radiation in the continuum does not depend on local 
properties (i.e., for an optically thin continuum, when it is given 
by the intensity emerging from the underlying photosphere). 
From the mathematical point of view, this method splits the 
Jacobi matrix into several parts and it is equivalent to a lin- 
earization of the whole system of model equations (i.e. hydro- 
dynamics, NLTE equations and radiative transfer). However, 
from the numerical point of view, it saves a huge amount of 
computer time because it is not necessary to invert large matri- 



ces. A similar method has been used, e.g., by Anderson J1985t . 
and Kubat (1994 1997) referred to it as to the implicit Un- 
earization of the ^j-factors. 

2.5. Calculation oftfie radiative force 

The radiative force is calculated using the Sobolev approxima- 
tion directly as the sum of contributions of individual lines after 



Castor ( fT974t . 



57r Vr \ ^ 



PiC^ r 



vH, j p{l + (jp^) (l-e-^-), (25) 



where is given by Eq. Ml\ . Occupation numbers Ni calcu- 
lated during the NLTE iteration step are used for a calculation 
of the radiative force. Note that in Eq. Ml\ number densities of 
individual levels are used instead of the total number density of 
a superlevel. These number densities may be calculated from 
Eq. ( I13> for those levels which are explicitly included in the 
NLTE rate equations. However, non-explicit levels, i.e., some 
of the higher excited levels or excited levels of the highest ion 
of an atom, are also included in the calculation of the radiative 
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force. For these non-explicit levels we use either the Boltzmann 
equilibrium at the wind temperature in the case of metastable 
levels or the radiative excitation equilibrium of Abbott & Lucy 
J1985> for other levels. The influence of these non-explicit lev- 
els on the radiative force is very small, however. 

The term (duri/dr) is calculated explicitly after 

i25\ where the variable a depends on a velocity gradient via 
i20\ . This term is included in the critical point condition (gen- 
eralised CAK condition, KKII, Eq. (57) therein). Derivatives 
of number densities calculated by means of Eq. (I24> are used 
to obtain derivatives of the radiative force with respect to 
fluid variables (densities, velocities, and temperatures) for the 
Newton-Raphson iterations. 

Oscillator strengths necessary for the calculation of the ra- 
diative force are extracted from the VALD database (Piskunov 
et al. 11995' Kupka et al. "1999V This set of data consists of 
about 180000 line transitions in the wavelength interval 250 - 
10 ODD A. 

2.6. Radiative cooling and tieating 

The radiative cooling and/or heating term is calculated using 
the thermal balance of electrons (Kubat et al. J_999 1. For the cal- 
culation of this term we include all the expUcit radiative bound- 
free and coUisional (both bound-bound and bound-free) transi- 
tions and free-free transitions of hydrogen and helium. Again, 
in order to obtain a well-converged model, derivatives of num- 
ber densities calculated by means ofEq. (Ell are used to obtain 
derivatives of the radiative cooling and/or heating terms with 
respect to fluid variables. The calculated derivatives of the ra- 
diative cooling/heating term are inserted into the Jacobi matrix 
of the Newton-Raphson iteration of the model hydrodynamic 
structure. Detailed expressions for the derivatives of heating 
and cooling rates are presented in Appendixicl 

2. 7. Electron density 

The charge of the passive component is calculated from ion- 
ization fractions of both hydrogen and helium. This enables us 
to calculate the correct value of the electron density, which is 
important for the proper solution of NLTE rate equations. 

2.8. Iteration scheme for the model calculation 

The initial estimate of fluid variables for the Newton-Raphson 
iterations (namely the so-called beta-law for velocities and den- 
sity obtained from the continuity equation) was described by 
KKII. For the initial value of the intensity we use an optically 
thin approximation, i.e., the radiation field is defined by the 
lower boundary condition at the beginning of our calculation. 

The computational procedure involves two similar basic 
steps, namely solution below and above the wind critical point. 
This splitting into two steps is necessary and is related to the 
critical point condition. In the first step the value of the base 
density, which allows the model to smoothly pass the critical 
point, is obtained. Afterwards, the model above this critical 
point is calculated. 



The calculation of a given model below the critical point 
consists of three nested iteration cycles (see Fig.Q. In the in- 
nermost iteration cycle the NLTE rate equations ([0 together 
with the radiative transfer equation both in lines (I14t and in 
continuum are solved. For this innermost iteration cycle the 
fluid variables (i.e. densities, velocities and temperatures of in- 
dividual wind components) are kept fixed. 

In the middle cycle we iterate the hydrodynamical structure 
of the model for fixed level populations. This iteration cycle is 
essentially the same as the iteration cycle of three-component 
nonisothermal wind models described by KKII. To save com- 
puting time, for large relative changes of fluid variables (greater 
than 10^^) we perform only 3 iterations of NLTE rate equa- 
tions, but when the relative changes of fluid variables are lower, 
and we may expect that we are closer to the correct solution, 
we perform more iterations of NLTE rate equations (up to 50). 
We assume that model is well converged if the relative changes 
of both fluid variables and occupation numbers are lower than 
10-5. 

In the outermost iteration cycle we solve for the value of 
the base density for which a model smoothly passes through 
the critical point. Again, these iterations were already described 
by KKII. For the base density (which determines the value of 
mass-loss rate dTl/dt) an accuracy of 1% is sufficient. 

After the base wind density is known, we calculate wind 
model above the critical point by a sequence of iterations of 
fluid variables. Each of these iterations is preceded by the set 
of NLTE iterations. Boundary conditions of the model above 
the critical point are specified at the critical point and are taken 
from the model below the critical point. The model below the 
critical point is fixed during the solution above the critical 
point. 

Splitting of a wind solution into two parts may cause 
problems since the radiation emitted from the outer region 
(i.e. above the CAK point) may influence the inner region. 
However, test calculations where we changed the outer bound- 
ary of inner solution showed that this is not the case, i.e. that 
the inner wind solution depends only marginally on the outer 
solution. The reason for such a behaviour is probably the strong 
dependence of the wind continuum optical depth on frequency. 
Calculations show that the stellar wind above the CAK point 
is essentially optically thin for continuum radiation with fre- 
quencies below the He II ionization limit, thus radiation is not 
absorbed in the continuum above the CAK point and does not 
influence the stellar wind below the CAK point. On the other 
hand, stellar wind above the CAK point is optically very thick 
for radiation with frequencies above the He II ionization limit. 
Again, due to this large optical thickness the radiation emitted 
in the regions above the CAK point does not influence the inner 
stellar wind. The influence of lines above the critical point on 
the subcritical region has been also neglected. 

3. Convergence behaviour of models 

As an example we study the convergence behaviour of a model 
corresponding to a Cam (HD 30614) with parameters given in 
Table|2l In Fig.|2lwe plotted relative change of population num- 
bers, and hydrodynamical variables (i.e. densities, velocities, 
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Fig. 1. Iteration scheme for the model calculation 



temperatures and charges of individual wind components) for 
the model below the critical point and for the model above the 
critical point. 

First we discuss the convergence behaviour of the model 
below the critical point. For the plot of relative change of level 
populations (upper left panel of Fig. ^ we combine NLTE 
iterations performed during individual hydrodynamical itera- 
tion steps for fixed base wind density (corresponding to the fi- 
nal CAK solution). Thus, in this plot the peaks corresponding 
to the hydrodynamical iteration steps occur They are caused 
by a reaction of NLTE populations to a significant change of 
the fluid variables after the calculation of the new estimate of 
the hydrodynamic structure. For the relative change of about 
10^^ relative changes oscillate. However, if we plot only rel- 
ative changes of populations of levels with largest occupation 
numbers (levels with Ni > 0.01), no oscillations appear. Thus, 



we conclude that the oscillatory behaviour is caused by weakly 
populated levels and that their influence on the convergence 
behaviour of the hydrodynamical structure is negligible (upper 
right panel of Fig.|3. This oscillatory behaviour is not a typical 
feature of convergence of our models. 

The plot of hydrodynamical structure convergence (upper 
right panel of Fig. |2} shows very fast convergence which is 
slightly slowed down by the oscillatory behaviour around the 
relative change 10^'^. This is caused by the change of the max- 
imum allowed number of NLTE iteration steps (see Sect. 12. 8> . 
Although due to this change the convergence is slightly slowed 
down, the total computing time is less because we do not lose 
time with calculation of precise occupation numbers for hy- 
drodynamical structure which is far from the converged one. 
Finally, note that the relatively small initial relative change of 
fluid variables (of order 10^^) of the hydrodynamical struc- 
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Fig. 2. The convergence of the models below (upper panels) and above (lower panels) the critical point. Upper left panel: The 
convergence of level populations below the critical point. Note that NLTE iterations in consecutive hydrodynamical iterations 
are combined in this plot. This explains the relatively large number of iterations and peaks that occur after the hydrodynamical 
iteration step. The oscillating behaviour for the relative change of order 10^^ is due to the weakly populated levels, because 
there are no oscillations in the plot of convergence of relative change of major levels (i.e. levels with Ni > 0.01). Upper right 
panel: The convergence of hydrodynamical structure (maximum relative change of all hydrodynamical variables, i.e. densities, 
velocities, temperatures and charges of individual wind components) of the model below the critical point. The slightly oscillating 
behaviour around the relative change 10^"^ is caused by the change of maximum allowed number of NLTE iteration steps (see 
the text). Lower left panel: The convergence of level populations of the model above the critical point. Again, NLTE iterations 
in different hydrodynamical iterations are put together and peaks occur after the hydrodynamical iteration step. Note the fast 
convergence of NLTE population between successive hydrodynamical iterations. Lower right panel: The convergence of the 
hydrodynamical structure of the model above the critical point. 



ture is caused by the way the critical point is calculated, ie. 
that models with similar base densities are calculated. Thus, 
the initial change of fluid variables of the final model below 
the critical point is very small. The convergence behaviour of 
previous models (i.e. models which do not correspond to the 
critical solution) is similar, however the initial change of fluid 
variables is higher. 

The convergence behaviour of the model above the critical 
point (lower panels of Fig.|2ji is similar to that below the criti- 
cal point. The convergence of occupation numbers (for a given 
hydrodynamical structure) is faster because the stellar wind is 
optically thinner (note that this fast convergence of NLTE iter- 
ations combined with larger number of iterations of hydrody- 



namic structure yields something like a "forest" of iterations in 
the plot, which is caused by a large change of occupation num- 
bers after a new hydrodynamic structure was calculated. The 
convergence of a hydrodynamical structure (lower right panel 
of Fig.O is however slower than below the critical point. This 
is partially caused by lower maximum change of wind param- 
eters allowed during Newton-Raphson iteration step. This con- 
straint of the maximum change of wind parameters is necessary 
for the stability of the solution (see Krticka|2003 1. 

The convergence of the model temperature is very fast in 
the case of a model below the critical point. The relative change 
of temperature is typically ten times lower than the relative 
change of the rest of the hydrodynamical variables in this case. 
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Fig. 3. Comparison of calculated and observed (taken from 
LSL, BG) terminal velocities. 
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Fig. 4. Comparison of calculated (crosses) and observed (taken 
from LSL, BG, errorbars) ratio of the terminal velocity to the 
escape velocity Vao/vesc- 



On the other hand, the convergence behaviour of a model above 
the critical point is dominated by the temperatures since the rel- 
ative changes of temperature are the highest. 



4. Comparison of calculated and observed 
quantities 

Because in future we intend to calculate mainly wind models 
of B stars, to test the reliability of our assumptions we calcu- 
lated NLTE wind models of stars with spectral type relatively 
close to B, i.e. later O stars. For this purpose we selected O 
stars with a spectral type 06 or later from the set of Lamers 
et al. 



4. 1. Terminal velocities 

Wind terminal velocities of galactic O stars can be measured 
with a typical uncertainty of 100 km s^^. Thus, in principle, 
they may provide a relatively strict test for any wind model. 
Unfortunately, because terminal velocities depend on surface 
escape velocities (see Castor et al. 1975.) the precise prediction 
of wind terminal velocity relies on the precise knowledge of 
stellar mass and radius. Consequently, the long-standing prob- 
lem of "mass discrepancy" in O stars (cf. Herrero et al. 1992! 
Lanz et al. I1996> complicates a straightforward comparison of 
calculated and observed terminal velocities. This is especially 
the case of stars HD 30614, HD 188209, and HD 218915, for 
which we used the same stellar parameters, but which have dif- 
ferent observed terminal velocities. 

The comparison of observed (taken from LSL and Bianchi 
& Garcia |2002 hereafter BG) and calculated terminal veloci- 
ties for model stars shown in Fig.|3]implies that our NLTE wind 
models are able to predict relatively correct values of wind ter- 
minal velocities. The fit is improved compared to our previ- 
ous calculations presented in KKII (see Table|2j. Note that our 
calculations consistently solved the problem of high theoret- 
ical terminal velocities obtained by LSL using the so called 
"cooking formula" of Kudritzki et al. ( 1989 1, see also Table 2 
and Fig. 12 in KKIL However, since our models include several 
simplifying assumptions, it is necessary to test this conclusion 
against more advanced models. 

It remains unclear why LSL obtained too high values of 
terminal velocities when theoretical calculations (KKII) with 
the same force multipliers yielded terminal velocities consis- 
tent with observations. Note that Pauldrach et al. (200 1'i also 
predicted correct terminal velocities, however for a different 
set of stars. 

The calculated terminal velocities are, in fact, not too dif- 
ferent from those calculated by KKII, although the terminal 
velocities presented by KKII were obtained simply using the 
force multipliers of Abbott (1982i. This relatively surprising 
conclusion can be attributed to the fact that the a parameter cal- 
culated properly using line statistics (Puis et al. 2000 1 is very 



II995I hereafter LSL). We included only stars for which 
reliabile determination of their mass-loss rate exists in the lit- 
erature. Stellar parameters (masses, radii and effective temper- 
atures) taken from LSL are listed in Table|2l 



close to the values calculated by Abbott J1982> . And because 
the terminal velocities depend on the a value, terminal veloci- 
ties calculated in this paper and those using Abbott's multipli- 
ers (KKII) do not differ significantly. However, properly cal- 
culated k parameters (Puis et al. '2000^ are significantly lower 
than those calculated by Abbott ( 1982 1. Indeed, our mass-loss 
rates (which are, in fact, influenced by the k parameter) are 
systematically lower than those of KKII and are close to the 
observed values (see also the next section for the discussion of 
mass-loss rates). 

The stellar wind theory predicts a strong correlation be- 
tween wind terminal velocity and stellar escape velocity. Thus, 
we have also compared the ratio of the observed and calculated 
terminal velocity to the escape velocity Uoo/wcsc (see Fig. |3 
Wcsc is calculated using stellar parameters from Table |3- Our 
derived mean ratio Woo/^csc ~ 2.6 is in a good agreement with 
the mean observed Voo/wcsc ratio 2.5 derived by LSL for stars 
with an effective temperature higher than 25 000 K. 
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Table 2. Stellar and wind parameters of selected O stars from our test sample. Stellar parameters (columns 971, i?*. Tog) were 
adopted from Lamers et al. ( 19951. Observed mass loss rates were taken from [1] - Leitherer ( 1988 ), [2] - Scuderi et al. J1992> . 
[3] - Lamers & Leitherer {19931, [4] - Puis et al. ( 1996 1, [5] - Crowther & Bohannan ( 1997 1, [6] - Scuderi et al. ( 19W), [7] 
- Lamers et al. 119991 [8] - BG (see the text), predicted values of Tl and v^o are determined using our code. Note that terminal 
velocities calculated by KKIl were obtained using simplified force multipliers after Abbott ( 1982). 
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Fig. 5. Comparison of calculated and observed mass-loss rates. Fig. 6. Comparison of theoretical mass-loss rates. 



4.2. Mass loss rates 

The credibility of observed mass-loss rates in most stars is 
much worse than their wind terminal velocities. Although in 
general the wind models are able to predict coiTect mass-loss 
rates (cf. Vink et al. '2000), there is a large scatter for indi- 
vidual stars. This situation is demonstrated in Fig. [S] where 
we compare calculated and observed mass-loss rates. We com- 
piled observed mass-loss rates based on the Ha measurement 
(Leitherer . 1988. Scuderi et al. [T992l Puis et al.[l996, Crowther 
& Bohannan IT9971 Scuderi et al. [T998l Lame rs et a l. [T999b 
and on the radio emission (Lamers & Leitherer [1993] Scuderi 
et al. I1998J together with mass-loss determinations obtained 
using UV lines by BG. Generally, UV line-diagnostic is be- 
lieved to be badly affected by the unknown wind ionization 
fractions (Lamers et al. il999j . This is especially the case of the 



O star winds for which unsaturated wind lines originate only 
from trace elements. However, precise modelling can overcome 
these difficulties (BG). 

We conclude that there is generally a good agreement be- 
tween predicted and observed mass-loss rates. However, sig- 
nificant differencies exist for individual stars. Similarly, com- 
parison of mass-loss rates calculated in this paper and by Vink 
et al. i2000> shows relative good agreement between these pre- 
dicted mass-loss rates. This possibly supports the conclusion 
of Puis ( 1987.) that multiple scattering (included into Vink et 
al. J2000> models) does not significantly influence mass-loss 
rates of weaker stellar winds. 

An important caveat has to be mentioned. The observed 
mass-loss rates may be influenced by wind clumping (e.g. 
Crowther et al. (TMT), Massa et al. (^UHT), Bouret et al. $^M^ ) 
and thus, the observed mass-loss rates shall be lower in such a 
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Fig. 7. Comparison of calculated (crosses) and observed (error- 
bars) wind momentum (determined by Herrero et al. QQQ2 1 for 
Cyg OB2 association). Linear regression of calculated data is 
also plotted. 

Table 3. Comparison of wind momentum-luminosity relation- 
ship from different sources, both theoretical and observational 
(see Eq. (I26» . Parameters of the sample denoted as "super- 
giants" were obtained by Repolust et al. (2003 1 combining their 
own results and results by Herrero et al. (2002i for supergiants. 
Similarly values of the sample "giants and dwarfs" were ob- 
tained. 



Sample 


log Do 


X 


supergiants 


17.98 ± 1.88 


2.00 ±0.32 


giants & dwarfs 


18.70 ± 1.29 


1.84 ±0.23 


Vink et al. ( 2000 1 


18.68 ±0.26 


1.826 ± 0.044 


this paper 


18.7 ±2.3 


1.83 ±0.40 



However, recent observational studies (Repolust et al. 
120031 Markova et al. I2003> indicate that the situation is prob- 
ably more complicated. They concluded that there are two 
different wind momentum-luminosity relationships for super- 
giants (or, probably more correctly, for stars with the Ha 
line in emission), and for other stars (mainly for giants and 
main-sequence stars). However, theoretical studies (cf. Vink et 
al. I2000> do not predict such a difference. Also results from 
our sample do not exhibit any significant difference between 
wind momentum-luminosity relationships for supergiant and 



for non-supergiant stars. Repolust et al. J2003> conclude that 
the observed difference may be a consequence of clumping and 
that observed mass-loss rates of high-density winds (mostly su- 
pergiant winds) may be artificially enhanced by clumping. To 
test this idea we calculated our predicted wind momentum- 
luminosity relationship and compared it with both the the- 
oretical one (Vink et al. I2000> and with the observed wind 
momentum-luminosity relationship for supergiants and non- 
supergiant stars (see Table|3}- Apparently, our predicted regres- 
sion to the linear wind momentum-luminosity relationship 

log [mv^ (i?/i?0)'/'] = xlog{L/LQ) + logi^o, (COS) 

(26) 

where x and Dq are fit parameters, agrees well with the theo- 
retical results of Vink et al. (2000 1 and with the observed rela- 
tionship for non-supergiant stars. However, the observed wind 
momentum-luminosity relationship for supergiants is different. 
This confirms the finding of Repolust et al. (2003j that the the- 
oretical wind momentum-luminosity relationship corresponds 
to the observed one for non-supergiant stars. Note, however, 
that most supergiants observed by Herrero et al. (.2002 ) repre- 
sent exceptions to this rule since most of them do not have an 
Ha line in emission. 



case (see also Kramer et al. 2003 for the discussion of effect of 
clumping on the theoretical profiles of X-ray lines). However, 
in such a case clumping would influence also the theoretical 
mass-loss rates (cf. Grafener l2003> . 

4.3. Wind momentum luminosity relationsiiip 

The wind momentum-luminosity relationship (see Kudritzki & 
Puis 1200 1 may provide a possibility to obtain precise mea- 
surements of stellar and galaxy distances. Thus, the exact cal- 
ibration of this relationship is necessary. Moreover, the wind 
momentum depends on the stellar mass only marginally, and 
thus it provides a more reliable tool for the investigation of the 
consistency between theoretical models and reality (Puis et al. 
[1996 1. 

We have compared a theoretical modified wind 

1/2 

momentum-luminosity relationship OJlwoo (R/Rq) 
obtained for considered stars with our NLTE wind models 
with a relationship derived from observations of supergiants by 
Herrero et al. ( I2002> for the Cyg OB2 association (see Fig.0. 
Again, there is a relatively good agreement between calculated 
and observed data. 



5. Discussion 

The main improvement of the presented models is that the 
radiative force is calculated without simplifying assumptions 
about the line-distribution function (Puis et al. 2000 1, i.e., the 
parameters k, a, and 5 (Castor et al. ll975l Abbott J_982,i are no 
longer required. A correctly-calculated wind temperature struc- 
ture influences occupation numbers and thus, also the radiative 
force. Moreover, these models are, in principle, able to consis- 
tently calculate a multicomponent wind structure, which is not 
necessary for the case of O star winds, but that will be invalu- 
able for calculations of low-density B star winds. 

On the other hand, there are several assumptions involved 
in the calculation of these models which can affect their valid- 
ity. 

First, we assume a stationary and spherically symmetric 
stellar wind. Both of these assumptions put severe limitations 
on our models. Variable phenomena, which are often observed 
in spectra of O stars, cannot be properly described with a sta- 
tionary spherically symmetric model of a wind. Spherical sym- 
metry also prevents a proper description of rotating winds. On 
the other hand, a detailed numerical study of wind stability (cf. 
Feldmeier at al. 119971 Runacres & Owocki '2002), and of ef- 
fects of stellar rotation on the stellar wind (cf. Petrenz & Puis 
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12000 1 showed that basic wind parameters (i.e. terminal veloc- 
ity and mass-loss rate) can be well reproduced by stationary 
and spherically symmetric wind models (for rotational veloci- 
ties lower than the critical rotational velocity). However, wind 
clumping may affect the values of theoretical mass-loss rates 
(Grafener l2003> . Finally, shocks caused by wind instabilities or 
by the presence of strong magnetic fields (ud-Doula & Owocki 
I2002> are sources of strong X-ray and EUV emission. This 
emission influences the ionization balance due to the Auger 
ionization (Cassinelli & Olson . 1979 1 and direct photoioniza- 
tion (Pauldrach et al. 1994). On the other hand, MacFarlane et 
al. (1994 1 showed, that for O stars the X-ray radiation influ- 
ences the ionization balance of trace elements only. Thus, we 
conclude that at least wind models of hotter star from our sam- 
ple are not significantly influenced by X-ray shock emission. 

Another sources of possible errors are included in the ap- 
proximate form of model equations itself. First, multiple scat- 
tering and line overlap may influence the mass-loss rate, espe- 
cially for high-density winds (see Abbott & Lucy 1985. Puis 
119871 Vink et al. [19991. However, we plan to apply our code 
first to models of low-density winds, for which the effect of 
multiple scattering is not important. 

Second, radiative transfer in lines may be influenced by the 
continuum effects. To overcome this limitation, we plan to in- 
clude continuum radiation in the formulation of the radiative 
transfer equation in the Sobolev approximation into our mod- 
els after Hummer & Rybicki ( 1985 1. For some lines, specially 
for the strong ones, the Sobolev approximation may not be 
credible. Thus, we plan to test the reliability of the application 
of the Sobolev approximation for the calculation of the radia- 
tive force by using a correct solution of the radiative transfer 
equation in moving media (Korcakova & Kubat |2C)03> . This 
will also allow as to account for the self-shadowing by photo- 
spheric lines which is important for thin winds (Babel . 1996t 
together with calculation of model atmospheres with inclusion 
of more elements). Moreover, it will enable a study of instabil- 
ities connected with the source-function gradients (Owocki & 
Puls [T999l . 

The effect of wind blanketing (Abbott & Hummer '1985 1 
was neglected in our calculations to make the treatment of the 
wind easier. However, it is clear that this effect causes back- 
warming of stellar surface and may lead to different values of 
effective temperatures (e.g. Herrero et al. 2002j. 

The lambda-iterations employed for the solution of rate 
equations together with radiative transfer in continuum may be 
source of another potentially very serious problem. However, 
the stellar wind is optically thin in most continuum frequencies 
and the optically thick region occurs mostly relatively near to 
the star. On the other hand, the radiative transfer in lines in a 
Sobolev approximation is solved consistently with rate equa- 
tions. Thus, due to these reasons and due to a relatively good 
convergence of our models we conclude that we are able to cal- 
culate consistent occupation numbers and wind structure even 
with this approximation. 

The processes of ionization and excitation can be very slow 
compared to the typical flow-time. Thus, the ionization state 
of the wind can be "frozen" in a gas parcel and, consequently, 
can be advected out. Such a case in the stationary spherically 



symetric stellar wind is described by the rate equation in the 
form of (Mihalas [T978t 



1 d (r'^riiVri) 



ni ^ Pi 



0, 



(27) 



instead of Eq.Q. However, numerical tests showed that this 
phenomenon can be neglected for most considered levels (see 
also Lamers & Morton 1976 1. Moreover, there would not be 
any radiatively driven wind if the opposite is true since radia- 
tively driven stellar wind is enabled by high transition rates Pij 
(Gay ley 1995 1. On the other hand, at first sight it is not clear 
whether this assumption is also fulfilled in the case of strong 
wind-shocks. 

Different approximations of NLTE model atoms may influ- 
ence the wind parameters obtained. Besides the well-known er- 
rors in the atomic databases, inclusion of superlevels may badly 
deteriorate calculated models. Whereas the effect of our sim- 
plified inclusion of superlevels is found to be small, neglected 
coUisional transitions between individual levels forming the su- 
perlevel may affect wind temperature structure (cf. Hillier & 
Miller [T998t . Moreover, another subtle basic effect may influ- 
ence the calculated occupation numbers. The calculated ion- 
ization fractions depend on the partition function approxima- 
tion. Especially, in the relatively thin wind environment (where 
a large number of higher excited levels have high occupation 
probabilities) different partition function approximations result 
in different final ionization fractions. Thus, it would be con- 
venient to include more sophisticated methods for the calcula- 
tion of partition functions (cf. Hummer & Mihalas[l988j in the 
model equations than those mentioned in Sect. 12.1. 11 

6. Conclusions 

We have presented new models of line-driven winds of hot 
stars. The radiative force and the radiative heating term are 
consistently calculated using statistical equilibrium equations. 
The statistical equilibrium equations are solved for 15 most 
abundant elements. Obtained occupation numbers are used for 
the calculation of the radiative force (in total about 180000 
lines are accounted for in the calculation of the radiative force) 
and for the calculation of radiative cooling and heating term. 
All bound-bound collisional transitions and bound-free coUi- 
sional and radiative transitions explicitly included in the mod- 
els are allowed to contribute to the cooling and heating terms. 
Moreover, free-free transitions of H and He are also included. 

We test the ability of our new models to predict correct 
wind parameters of late O stars. In particular, we compare 
calculated terminal velocities, mass-loss rates and wind mo- 
menta with the observed ones. We show that wind parameters 
calculated using our models are consistent with observation. 
Moreover, our new wind models have solved the problem in- 
troduced by LSL of high theoretical wind terminal velocities. 
However, this conclusion has to be tested again against more 
advanced models, because our models involve several simpli- 
fying assumptions. Moreover, the relatively good correspon- 
dence between observed and theoretical wind parameters can 
be reduced by our insufficient knowledge of stellar and wind 
parameters. 
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There is a good agreement between our calculated wind 
momentum-luminosity relationship and that relationship ob- 
served by Herrero et al. l l20(32t . and Repolust et al. (2003 1 for 
non-supergiant stars. However, contrary to the theoretical pre- 
diction of Vink et al. 12000) and ours, there is a difference be- 
tween the wind momentum-luminosity relationship observed 
for non-supergiant stars and for supergiants (or, more precisely, 
between stars with a different type of Ha profile). Repolust et 
al. (|2003 1 attribute this difference to wind clumping. 

On the other hand, many steps remains to be done. Besides 
the overall improvements of our code we have to compare ob- 
served and calculated wind parameters for B stars before we 
start to study multicomponent models. 

We postpone detailed discussion of wind temperature struc- 
ture for a separate paper 
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Appendix A: Atomic data for NLTE calculations 

The source of all data used for the calculation of model atoms is 
given in Table lA.ll Most important models are taken from the 
TLUSTY web page |http : //tlusty . qsf c . nasa . qov| 

These atomic models were complemented by data con- 
structed for some light elements using the Opacity Project 
(hereafter OP; Seaton [T987l Luo & Pradhan 1989 Seaton et 
al. [T992l Butler et al. [T993l Nahar & Pradhan 1993 1 database 
fhttp :/ / vizier . u-strasbg . f r/OP . html Energy 
levels, photoionization cross sections and oscillator strengths 
were taken from OP data. To avoid complicated calculation 
of numerous resonances (which shall be, in fact, accounted 
for using solution of radiative transfer equation in moving 
medium), we applied linear regression of photoionization 
data. For the collisional ionization we applied the so-called 
Seaton's formula (Mihalas 11978] Eq. (5.78) therein) with the 
photoionization cross section at the ionization edge taken from 
OP data. Finally, for the collisional excitation we used the Van 
Regemorter formula (Mihalas 11978! Eq. (5.75) therein. Van 
Regemorter 19601. 

Another important atom is iron. For Fe III we used energy 
levels and oscillator strengths from Nahar & Pradhan (1996), 
photoionization cross sections from OP (Sawey & BeiTington 
|1992) , collision rate coefficients for low-lying levels from 



Zhang ( I1996> . and the Van Regemorter formula for levels not 
included in Zhang (19961 Radiative rates for Fe IV and Fe V 
were calculated using data from Bautista & Pradhan Cl997j or 
Bautista 



I1996> . and collision excitation rates are again calcu- 
lated using either data from Zhang & Pradhan (1997) or the 
Van Regemorter formula. Finally, radiative cross sections (both 
bound-bound and bound-free) for Fe VI are taken from OP data 
and collisional rate coefficients are taken from Chen & Pradhan 



Appendix B: Details of linearization of rates 

Solution of NLTE rate equations is performed using Newton- 
Raphson iterations (the so called linearization). Linearization 
matrix follows from Eq. (I22> . where derivatives of corre- 
sponding to the line transitions are evaluated using Eqs. (|3}, 
(13, and (El) - (EU as 
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Note that the derivatives of collisional rates with respect to Ni 
are zero. Finally, the derivatives of normalization condition (|5J 
with respect of all relative populations Ni are 
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dN 



(B.3) 



where k is index of level with the highest occupation number. 

Similarly, for the calculation of derivatives of relative num- 
ber densities N with respect to fluid variables Uc, Pp, To, cr and 
Vri in Eq. (I24t it is necessary to evaluate quantities dPij/dh. 

The derivatives of excitation and deexcitation rates with re- 
spect to electron density Ue follow from Eqs. (0, (EJ, and (|8}, 
and are relatively simple. 
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Table A.l. Details of atomic models for NLTE calculation 
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whereas for the calculation of the derivatives of the ionization 
rates (Eqs. and (|8}) it is necessary to account for the depen- 
dency of Saha-Boltzmann ratio on the electron density 
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The derivatives of the excitation rates with respect to elec- 
tron temperature are also relatively simple. 
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however, for the evaluation of derivatives of ionization rates it 
is necessary to perform an integration 
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Derivatives of collisional rates are given by Eqs. ( IB.6bt and 

Using the dependence of the radiative transfer equation in 
Sobolev approximation on the velocity gradient a and on Vri 
and pp via the fraction y = Pp/vri we calculate derivatives of 
excitation rates with respect to these variables, namely 
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Note that the variable cr is not an indepentend variable during 
linearization. Derivatives with respect to a are transformed to 
derivatives of ionic velocity using Eq. MQ\ . 

Appendix C: Details of iinearization of raditive 
heating/cooiing term 

Radiative heating/cooling term in the case of thermal balance 
of electrons is given by (Kubat et al. ll999> 



where individual terms are 
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where ki denotes level into which is level i ionized. From these 
equations derivatives of heating/cooling term with respect to 
Vri, Pp, Pc,Tc for the Newton-Raphson iteration step of hydro- 
dynamical structure can be calculated. 
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